---
title: "Tables"
author: "Lukas Fesenfeld"
date: "2/12/2021"
output: html_document
---
```{r, echo=FALSE, message = F}
#run this the first time you run the code:
    #devtools::install_github('junkka/ehahelper')

library(here)
library("ehahelper")
library("broom")
library("knitr")
library("kableExtra")
library(tidyverse)

load('code_12_02_2021.RData')

options(digits = 3)

```

### Body Tables
```{r, echo=FALSE, message = F}

  t1 <- tidy(cox.mod3a_control_full_rand, exp = T) #convert regression output to dataframe, setting exp argument to "true" gives hazard ratio
  t1_clipped <- t1[c(1:2),] #only show relevant variables
  colnames(t1_clipped)[1] <- " "
  colnames(t1_clipped)[2] <- "hazard ratio"
  colnames(t1_clipped)[4] <- "z-statistic"
  
  #create table
 table1 <-  t1_clipped %>%
    kbl(caption = "<center><strong>H1a & H1b: Using coxme random intercepts per country</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")
 
 table1
```

```{r, echo=FALSE, message = F}
#h2a & h2b: using coxme random intercepts per country
  
  t2 <- tidy(cox.mod3_control_full_rand_high, exp = T)#convert regression output to dataframe
  t2_clipped <- t2[c(1:2),] #only show relevant variables
  colnames(t2_clipped)[1] <- " "
  colnames(t2_clipped)[2] <- "hazard ratio"
  colnames(t2_clipped)[4] <- "z-statistic"
  
  #create table
  table2 <- t2_clipped %>%
    kbl(caption = "<center><strong>H2a & H2b: Using coxme random intercepts per country</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")
  
  table2
```

```{r, echo=FALSE, message = F}
#h3a & h3b: using coxme random intercepts per country
  
  t3 <- tidy(cox.mod3_control_full_rand_low, exp = T)#convert regression output to dataframe
  t3_clipped <- t3[c(1:2),] #only show relevant variables
  colnames(t3_clipped)[1] <- " "
  colnames(t3_clipped)[2] <- "hazard ratio"
  colnames(t3_clipped)[4] <- "z-statistic"
  
  #create table
  table3 <- t3_clipped %>%
    kbl(caption = "<center><strong>H3a & H3b: Using coxme random intercepts per country</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")
  
  table3
```

### Visualisation
```{r, echo=FALSE, message = F, fig.height = 5, fig.width = 7, fig.align = "center"}
#plots

# The first figure would only show the effects of IPA on the three outcome variables of interest, i.e. the hazard rate of experiencing a) discontinuation, b) supersession with a policy with a lower IPA value, c) supersession with a policy with a higher IPA value



plotdata <- rbind(t1_clipped, t2_clipped, t3_clipped)
plotdata <- plotdata[c(1,3,5),c(1,2,6,7)]
plotdata$outcome <- factor(c("Discontinuation",
                      "Supersession with a policy \n with a higher IPA value",
                      "Supersession with a policy \n with a lower IPA value"), 
                      levels = c("Supersession with a policy \n with a lower IPA value",
                                 "Supersession with a policy \n with a higher IPA value",
                                 "Discontinuation"))

plotdata2 <- rbind(t1_clipped, t2_clipped, t3_clipped)
plotdata2 <- plotdata2[c(2,4,6),c(1,2,6,7)]
plotdata2$outcome <- factor(c("Discontinuation",
                      "Supersession with a policy \n with a higher IPA value",
                      "Supersession with a policy \n with a lower IPA value"), 
                      levels = c("Supersession with a policy \n with a lower IPA value",
                                 "Supersession with a policy \n with a higher IPA value",
                                 "Discontinuation"))
                      



plot1 <- ggplot(data = plotdata, aes(x = `hazard ratio`, y = outcome))+
  geom_point()+
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, height = 0))+
  labs(title = "Effects of IPA on outcome variables of interest (Hazard Ratio)",
       subtitle = "95% confidence intervals",
       x = "Estimate",
       y = "")+
  xlim(0,5)+
  geom_vline(xintercept = 1, linetype = "dotted")
  

plot2 <- ggplot(data = plotdata2, aes(x = `hazard ratio`, y = outcome))+
  geom_point()+
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, height = 0))+
  labs(title = "Effects of Technology Specificity on\noutcome variables of interest (Hazard Ratio)",
       subtitle = "95% confidence intervals",
       x = "Estimate",
       y = "")+
  xlim(0,2)+
  geom_vline(xintercept = 1, linetype = "dotted")


plot1
plot2

```


### Appendix Tables
```{r, echo=FALSE, message = F, results = F}
#get confidence intervals for summaries
fc <- c(NA, NA)

    cis11 <- as_tibble(summary(cox.mod3a_control_fixed)$conf.int[,3:4])
      colnames(cis11) <- c("conf.low", "conf.high")
     # cis11 <- rbind(cis11, fc)
      
    cis12 <- as_tibble(summary(cox.mod3a_control_full)$conf.int[,3:4])
      colnames(cis12) <- c("conf.low", "conf.high")
      cis12 <- rbind(cis12, fc)
    
    cis21 <- as_tibble(summary(cox.mod3_control_fixed_high)$conf.int[,3:4])
      colnames(cis21) <- c("conf.low", "conf.high")
     # cis21 <- rbind(cis21, fc)
    
    cis22 <- as_tibble(summary(cox.mod3_control_full_high)$conf.int[,3:4])
      colnames(cis22) <- c("conf.low", "conf.high")
      cis22 <- rbind(cis22, fc)
    
    cis31 <- as_tibble(summary(cox.mod3_control_fixed_low)$conf.int[,3:4])
      colnames(cis31) <- c("conf.low", "conf.high")
    #  cis31 <- rbind(cis31, fc)
      
    cis32 <- as_tibble(summary(cox.mod3_control_full_low)$conf.int[,3:4])
      colnames(cis32) <- c("conf.low", "conf.high")
      cis32 <- rbind(cis32, fc)
      
    

#Model summaries
    sum11 <- tidy(cox.mod3a_control_fixed)  #fixed effects
    haz11 <- tidy(cox.mod3a_control_fixed, exp = T) #hazard ratio
    haz11 <- haz11[2]
    sum11 <- cbind(sum11[1:2], haz11, sum11[3:5], cis11 )
    colnames(sum11) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
   
    
    sum12 <- tidy(cox.mod3a_control_full)  #frailty
    haz12 <- tidy(cox.mod3a_control_full, exp = T) #hazard ratio
    haz12 <- haz12[2] 
    sum12 <- cbind(sum12[1:2], haz12, sum12[3:5], cis12)
    colnames(sum12) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    
    sum13 <- tidy(cox.mod3a_control_full_rand) #coxme random intercepts per country
    haz13 <- tidy(cox.mod3a_control_full_rand, exp = T) #hazard ratio
    haz13 <- haz13[2]
    sum13 <- cbind(sum13[1:2], haz13, sum13[3:7])
    colnames(sum13) <- c("", "estimate", "hazard ratio", "std.error", "z-statistic", "p.value","conf.low", "conf.high")
   
    sum21 <- tidy(cox.mod3_control_fixed_high)  #fixed effects
    haz21 <- tidy(cox.mod3_control_fixed_high, exp = T) #hazard ratio
    haz21 <- haz21[2]
    sum21 <- cbind(sum21[1:2], haz21, sum21[3:5], cis21 )
    colnames(sum21) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    sum22 <- tidy(cox.mod3_control_full_high)  #frailty
    haz22 <- tidy(cox.mod3_control_full_high, exp = T) #hazard ratio
    haz22 <- haz22[2]
    sum22 <- cbind(sum22[1:2], haz22, sum22[3:5], cis22 )
    colnames(sum22) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    sum23 <- tidy(cox.mod3_control_full_rand_high)  #coxme
    haz23 <- tidy(cox.mod3_control_full_rand_high, exp = T) #hazard ratio
    haz23 <- haz23[2]
    sum23 <- cbind(sum23[1:2], haz23, sum23[3:7])
    colnames(sum23) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    sum31 <- tidy(cox.mod3_control_fixed_low)  #fixed effects
    haz31 <- tidy(cox.mod3_control_fixed_low, exp = T) #hazard ratio
    haz31 <- haz31[2]
    sum31 <- cbind(sum31[1:2], haz31, sum31[3:5], cis31 )
    colnames(sum31) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    sum32 <- tidy(cox.mod3_control_full_low)  #fixed effects
    haz32 <- tidy(cox.mod3_control_full_low, exp = T) #hazard ratio
    haz32 <- haz32[2]
    sum32 <- cbind(sum32[1:2], haz32, sum32[3:5], cis32 )
    colnames(sum32) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")
    
    sum33 <- tidy(cox.mod3_control_full_rand_low)  #coxme
    haz33 <- tidy(cox.mod3_control_full_rand_low, exp = T) #hazard ratio
    haz33 <- haz33[2]
    sum33 <- cbind(sum33[1:2], haz33, sum33[3:7])
    colnames(sum33) <- c("", "estimate","hazard ratio", "std.error", "z-statistic", "p.value", "conf.low", "conf.high")


     


```

```{r, echo=FALSE, message = F}
#Create appendix tables: main models of interest

ms11 <- sum11 %>% kbl(caption = "<center><strong>Outcome variable: Discontinuation of policies <br>
Type of analysis: Cox proportional hazard analysis with country fixed effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms12 <- sum12 %>% kbl(caption = "<center><strong>Outcome variable: Discontinuation of policies <br>
Type of analysis: Cox proportional hazard analysis with frailty and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms13 <- sum13 %>%
  kbl(caption = "<center><strong>Outcome variable: Discontinuation of policies <br>
Type of analysis: Cox proportional hazard analysis with ‘coxme’ random effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")


ms21 <- sum21 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a higher IPA value (i.e., ratcheting-up) <br>
Type of analysis: Cox proportional hazard analysis with country fixed effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms22 <- sum22 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a higher IPA value (i.e., ratcheting-up) <br>
Type of analysis: Cox proportional hazard analysis with frailty and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms23 <- sum23 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a higher IPA value (i.e., ratcheting-up) <br>
Type of analysis: Cox proportional hazard analysis with ‘coxme’ random effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")


ms31 <- sum31 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a lower IPA value (i.e., ratcheting-down) <br>
Type of analysis: Cox proportional hazard analysis with country fixed effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms32 <- sum32 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a lower IPA value (i.e., ratcheting-down) <br>
Type of analysis: Cox proportional hazard analysis with frailty and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

ms33 <- sum33 %>% kbl(caption = "<center><strong>Outcome variable: Supersession with a policy with a lower IPA value (i.e., ratcheting-down) <br>
Type of analysis: Cox proportional hazard analysis with ‘coxme’ random effects and covariate-by-time interactions</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")
```

```{r, echo=FALSE, message = F}
#create appendix tables: robustness checks
zph1 <- readRDS("zph1.rds")
zph2 <- readRDS("zph2.rds")
zph3 <- readRDS("zph3.rds")

rt1 <- as.data.frame(zph1$table) %>% kbl(caption = "<center><strong>Robustness checks: H1</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

rt2 <- as.data.frame(zph2$table) %>% kbl(caption = "<center><strong>Robustness checks: H2</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")

rt3 <- as.data.frame(zph3$table) %>% kbl(caption = "<center><strong>Robustness checks: H3</strong></center>") %>%
    kable_classic(full_width = F, html_font = "Cambria")
```


```{r, echo=FALSE, message = F}
#print tables
ms11
ms12
ms13

ms21
ms22
ms23

ms31
ms32
ms33

rt1
rt2
rt3

```

